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ABSTRACT 

A technique is described for the detection and measurement of close binary 
systems whose images are unresolved. The method is based on analysis of the 
moment of inertia tensor of the image, from which the product of the binary 
flux ratio and square of the angular separation may be determined. Intrinsic 
asymmetries of the point-spread function are removed by comparison with the 
image of a reference star. Multiple exposures may be used to increase the 
signal-to-noise ratio without need of image alignment. An example is given of a 
simulated measurement of the dwarf carbon star system G77-61. 



1. INTRODUCTION 

There are a number of astrophysical situations where one needs to measure either the 
position or the flux of a faint unresolved object which is located close to a bright star. 
Obvious examples are the search for extra-solar planets and low-luminosity companions to 
bright stars, and the study of close binary systems. If the flux ratio between the primary 
and secondary object is very large, coronographic techniques may be employed, but this 
is can only be effective if the separation of the images of the primary and secondary 
objects is substantially larger than the radius of the point-spread function (PSF). For small 
separations alternate techniques are needed. 

When the separation between the binary components is comparable to, or smaller 
than, the radius of the PSF, the image of the system will be elongated to some degree. This 
paper describes a technique which uses the image elongation to obtain information about 
the photometric and astrometric parameters of the system. 



2. THE ALGORITHM 



There are a several conditions that are desirable in any technique which aims to 
measure image elongation. Some important ones are: 
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1. invariance to the location of the PSF 

2. invariance to the orientation of the detector 

3. insensitivity to the shape of the PSF 

4. optimization of the signal-to-noise ratio 

5. ability to co-add images to improve the signal-to-noise ratio 

The method described in this paper satisfies, to a large degree, all of these conditions. 



2.1. Characterizing Image Elongation 

We begin by reviewing some elementary properties of the moments of the light 
distribution. For a Cartesian coordinate system x = (x^,a;^) and intensity distribution 
/(x), the zeroth, first and second moments are defined by 



p = j f{^)d''x (1) 
/ = y f{:ic)x'(fx (2) 



f{x)x'x^(fx (3) 



where i — 1,2 and the function /(x) is presumed to be vanishingly small at large values 
of |x|. The integrals extend over a sufficiently large range that contributions from the end 
points are negligible. 

A shift in the position of the image, /(x) — > /(x — a), induces the following changes, 
which are easily obtained by the substitution y = x — a, 

P ^ P (4) 
p' p' + a'p (5) 
p^j ^ + ay + + aVp (6) 

from which it follows that the inertia tensor 

= Ip'^ - py/p]/p (7) 



is invariant under translations of the image. 
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Under rotations of the coordinate system, x'"^ — A*^x^ (a summation over repeated 
indices is implied), it is evident from the definition that p and transform as a scalar and 
vector respectively, and that both p^^ and transform as second-rank tensors, eg. 

V'^ = A\A'il''^ . (8) 

Where there is no ambiguity, we may omit the indices on vectors, tensors and matrices and 
use bold face type to distinguish them from scalars quantities. 



2.2. Moments of a Binary Image 

Let /(x) be proportional to the PSF. We chose the proportionahty constant and the 
origin of the x coordinate system in such a way that the zeroth moment is unity and the 
first moments vanish, ie 

P = 1 (9) 
p' = (10) 

where the symbol ' " ' indicates that these are normalized moments of the PSF. The second 
moments are generally nonzero and contain information about the size and shape of the 
PSF. 

Now consider an image of the form 

/(x) = /(x + a/2) + 6/(x-a/2) (11) 

where b < 1. This corresponds to a binary system where the secondary has a flux b, relative 
to the primary, and a vector separation a. Using Eqns. 1-3 and 7, we obtain the inertia 
tensor of this image 

I- = + (Yl^a'a^- ■ (12) 

From this we see that the presence of the secondary component adds a term to the inertia 
tensor of the PSF which, to first order, is proportional to the relative flux times the square 
of the separation. By subtracting the components of the inertia tensor of the PSF from 
that of the binary system, we obtain this term, which we denote by M, 

From this, if the separation vector a is known, we obtain the relative flux. In fact, only the 
magnitude a of the separation vector is needed - by taking the trace of the matrix M we 
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obtain 

^'■(M) = . (14) 

from which b may be determined. 

Alternatively, if the relative flux is known, the diagonal terms of M give the components 
of separation vector. Even if b is not known, the position angle of the separation vector, 
with respect to the axis, can be found from the ratio of diagonal terms, 

tan0 = aVa^ = (M'Vm")^/^ _ ^^5^ 



2.3. Sampling Effects 

In practice we work with a digital image which consists of a discrete set of samples 
fa corresponding to pixels located at positions x^, normally on a square grid at intervals 
p. The values /„ are proportional to the product of the intensity and the pixel response 
functions gQ.(x) integrated over the pixel areas. If the response function is the same for all 
pixels we have 



fa = J f{^)qa{^)d'x (16) 

= y'/(x)g(x-x«)rf2a; . (17) 



Thus the values fa actually sample the image formed by smoothing the incident intensity 
with the pixel response function. This effect can be incorporated by employing the actual 
PSF of the sampled image, which includes the result of the smoothing. 



From the sampling theorem (see, for example, [Bracewell 1978[ ) the function /(x) can 



be uniquely reconstructed from the samples fa only if it contains no spatial frequencies 
higher than the Nyquist frequency 1/p. For a diffraction limited image, at wavelength A, 
from a telescope of aperture diameter D, the pixel spacing must be no greater than X/2D. 
The image is then given by 



/(x) = /„sinc(x^ - a;^)sinc(x^ - xl) 

a 

where sinc(x) = sin(7ra;)/7rx. Substitution of this equation into Eqns. 1-3 gives the 
corresponding discrete expressions 



P = I] /a , 

a 



(19) 
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P' = T.fc^<^ (20) 

a 

P'' = Jlfc^^a^ (21) 

a 

providing that the integrals converge. Eqns. 7, 13-15 and 19-21 then provide a direct path 
from the observational data to the parameters b and a of the binary system. 

While this simple technique seems plausible, in practice it does not work due to 
problems of convergence and noise, discussed next. However a modification of the technique, 
developed in Section 2.5 avoids these difficulties while at the same time optimizing the 
signal-to-noise ratio. 



2.4. Convergence and Noise 

The above algorithm fails because for real astronomical images the integrals in Eqns. 
1-3 do not converge. This is a result of the PSF of an optical system being a band-limited 
function. The amplitude of the light in the PSF is the two-dimensional spatial Fourier 
transform of the optical transfer function (OTF). There is an upper limit to the spatial 
frequency at which the OTF can be non-zero, imposed by the finite size of the telescope 
entrance pupil. The OTF is the product of the instrumental and atmospheric transfer 
functions with the pupil function, defined as unity for points within the pupil and zero for 
points outside. From the convolution theorem, the amplitude of the PSF is therefore the 
convolution of the instrumental and atmospheric response with the Fourier transform of the 
pupil function. Since the square of the latter function is just the diffraction-limited PSF, it 
follows that the intensity of light in the PSF cannot fall off faster than the diffraction limit. 
For a circular aperture this has the form (eg. Porn fc Wolfe 1980|) 



/(x) = [2Ji{kRr)/kRrf (22) 

where r = |x| is measured in radians, k = 2tt/X is the wave number, A the wavelength, and 
R = D/2 is the radius of the aperture. For large values of its argument, the Bessel function 
has the asymptotic expansion 

Ji(x) — > — (sinx — cosx) (23) 
so / falls off no faster than r^^. Since d'^x cx r, the integral in Eqn. 3 is ill-defined. 



Even more serious is the fact that random noise, present in the image, causes a rapid 
divergence of the moments. To see this, recall that both the photon noise and the read 
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noise is uncorrelated from pixel to pixel. The presence of such noise adds a fluctuating 
component to /(x) which in turn produces fluctuations in the moments. The variance of 
these fluctuations is the sum of the variances of the individual pixel fluctuations, weighted 
by the squares of the coordinate terms present in Eqns. 2 and 3. Since the read noise is 
independent of position, the variance of all moments increases without limit over the range 
of integration. 



2.5. Weighted Moments 

The solution that we propose for these problems is to apply weight factors to the 
intensity values of the image before computing the moments. In order to preserve 
the transformation properties of the inertia tensor, the weights cannot depend on the 
coordinates, but can only be a function of /„. We make the substitution — > Wafa, and 
ask what choice of Wa minimizes the random error in the moments. Prom Eqns. 19-21 we 
then have 

Var{p) = Y^Varifa) (24) 

a 

Var{p') = 5:(a;L)Var(/«) (25) 
Var{p^^) = Eix^xifVarif^) (26) 

a 

where Var{x) denotes the variance of the random variable x. We seek to minimize the 
ratios Var{p)/{p)'^, Var{p^)/{p^Y and Var{p^^)/{p^^y, and so equate to zero the derivatives 
of these ratios with respect to Wa- This leads to the result 

- (27) 



Two possible dependencies for the noise are Var{fa) oc and Var{fa) — constant. The 
first case corresponds to photon (Poisson) noise. It leads to = constant which, as we 
have seen, is unacceptable. The second case corresponds to a constant read noise. It leads 
to the choice Wa oc fa- As the inertia tensor is independent of the normalization of fa, 
we may choose the proportionality constant to be unity. This choice of weights solves the 
convergence problems as the asymptotic dependence of Wafa = fa ^^'^ proportional to 



While the use of weights solves the convergence and noise problems, it introduces 
complications in the interpretation of the results. Our choice of weights is equivalent to 
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performing the previous (unweighted) analysis on a new image formed by squaring the 
intensity values of the original image. Thus Eqn. 11 is replaced by 

/(x) = [/(x + a/2) + 6/(x - a/2)]2 . (28) 

After some algebra, we obtain the inertia tensor 

_ ^ b ( (1 + 6^)7(a) + 26 1 2bj{a) ^ _ ^ 

^ +2\[l + 62 + 267(a)P/"" +l + 267(a) + 62L^ ^^^^ 

where 

//(x + a/2)/(x-a/2)d2x ^ ^ 

7 a = '-^r^ '—^ 30 

//2(x)d% ^ ' 

pj(. ^ //(x + a/2)/(x-a/2)xVd^x 

This result is similar in form to Eqn. 12 but has an extra term involving the factor 
/'•'(a) — P^. This factor depends on the shape of the PSF and the separation a. For 
any separation, it can be computed numerically, using the image of the reference star to 
estimate the PSF. 

The function 7(a) approaches unity for small values of a and zero for large a. For a 
Gaussian PSF it has the simple analytic form 

7(a) = exp(-aV4(T2) (32) 

where a is related to the FWHM w of the PSF by w'^ — 8 In 2(7^. For close binary systems 
a/ a is small and 7(a) is quite close to unity. 

Eqns. 29-31, which represent the generalization of Eqn. 12 in the weighted case, relate 
the binary system parameters to the difference in the inertia tensor components between 
the binary and comparison images. 



2.6. Co-adding Images 

The signal-to-noise ratio that can be achieved in a single image is limited by the 
dynamic range of the detector. However, the signal-to-noise ratio can be increased by the 
use of multiple exposures. A high-signal-to-noise-ratio image can be obtained by co-aligning 
the individual images, interpolating and summing the individual intensities. However, the 
complexity of this process can be avoided if one wants only to measure the image ellipticity. 
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The moments, Eqns. 1-3, are linear in /. Because of this, the sum of the moments 
of a number of individual images are equivalent to the moments of the summed image. 
Moreover, because I is independent of translations, it is not necessary or desirable to 
co-align the individual images. One simply computes the inertia tensor Iq, for each image 
and averages the corresponding components: 

1 ^ .. 

= 77 E 4^ • (33) 

This method is preferable to computing the inertia tensor for a coadded image because the 
latter technique is sensitive to errors in alignment of the individual images. 



3. PERFORMANCE 

3.1. Noise Analysis 

The analysis of noise propagation in Equations (7), (19-21), and (28) is straight 
forward. Considerable simplification results if wc make the reasonable assumption that the 
PSF is circularly symmetric, at least to first order. For the limiting cases in which the noise 
variance per pixel is cither constant (read noise dominated) or proportional to the intensity 
(photon noise dominated) we find 

Varin ^ \[{n' + 1^] (34) 
where s is the total signal-to-noise ratio of the image, 

. = p/[Var{pp' (35) 

and 

jijkl _ Z^yJa) •^a-^a-^a-^a (op\ 

" - E(/a)"/^ ^ ^ 

where n = 2 for the read-noise-dominated case and n = 3 for the photon-noise-dominated 
case. 

For the case of diffraction-limited imaging by a circular aperture, the tensors may be 
evaluated analytically. Using Eqn. 22 and replacing the summation by integration in Eqn. 
36, and recalling that we are working with the square of the intensity (Section 2.5), we 
obtain 

.11^.22 ^ 16 /o^" S^[J,{kRr)/kRrYcos\4>)T^dTd4> 
327r j^[Ji{kRr) /kRrYrdr 
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1 POO poo 



2F 

= 0.088179(/ci?)-2 . 



(37) 



Similarly, 



7-1111 _ 7-2222 _ 7-1122 
-'2 — ^2 — ^2 



7-1111 _ 7-2222 _ 7-1122 
-'3 — -'3 — -'3 



2 roo roD 

Mxrxdx/j^ Mxrx-'dx 

9.76197(A;i?)-^ , 

-—- / JAxfx-^dxl / JAxfx-^dx 
8k^R^ Jo ^ ' 'Jo ^ ^ 

1.01967(A;r)-^ 



which gives 



Var{I^^) = 



10.53953 ' 




A 


1.797219 ^ 




27rR_ 



9.76197 ' 




A 


1.01967 ^ 







(38) 
(39) 

(40) 
(41) 



where the upper numerical values refer to the read-noise-dominated case and the lower 
values refer to the photon-noise-dominated case. 

Prom Eqn. 29 it is evident that there is no simple relation between the binary system 
parameters and the inertia tensor. Simulations indicate that the last term in this equation 
typically contributes about 1/3 of the ellipticity signal. If we assume that the variance in 
the image of the comparison star image is comparable to that of the binary system, we 
expect that 

Var{ba'a^) - 10 Var{P^) . (42) 

Eqns. 40-42 can be used to estimate the signal-to-noise ratio that would be required for 

any desired measurement of the binary system parameters. For example, if we assume that 
the magnitude a of the separation vector is known, the error in the measured flux ratio 
becomes, for the photon-noise-limited case. 



a{b) = [Varib)] 



1/2 



1 



A 

oD 



(43) 



Thus, in order to detect and measure, to 10% accuracy, a binary which has flux ratio 0.01 
and separation comparable to the size of the diffraction-limited PSF, a — X/D, we need 
a{b) — 0.001 which implies a total signal-to-noise ratio of order s ~ 10^. 
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3.2. Simulations 

In order to verify this analysis and evaluate the performance of the method, numerical 
simulations were performed. In a series of Monte-Carlo runs, artificial images were created 
of both a binary system and a reference star, by placing a suitably-scaled PSF image at 
the location of each component and then adding random noise. For each run the program 
estimated the desired parameters. The mean values and standard errors of the parameters 
were then determined for the set of runs. 



As an example, consider the dwarf carbon star system G77-61 ( Pahn et al. 197? 



Dearborn et al. 1986| ). The secondary is believed to be a population-II white dwarf, and a 



measure of its luminosity can provide a lower limit to the age of the universe ( [Richer et al. 



19971) . In the near infrared, the white dwarf is expected to be 4 to 5 magnitudes fainter 
than the primary star, and the maximum separation of the system is estimated to be 0.035 
arcsec. For a hypothetical observation using the Hubble Space Telescope (HST) and NICl 
infrared camera. The pixel size is 0.043 arcsec and the aperture diameter is 2.4 m. Thus, to 
avoid aliasing, we must use filters which block wavelengths below 1.0 um. The flux ratio is 
expected to be 0.009 and 0.021 at a wavelength of 1.10 um and 1.55 um respectively. We 
regard the separation as known and wish to measure the relative flux from the observations. 

The results of simulations of this system are presented in Figs. 1 and 2 and Table 1. 
These show the estimated value of the flux ratio b, and its statistical error, as a function 
of the total signal-to-noise ratio s, for two different wavelengths. At low values of s, noise 
fluctuations dominate, producing a substantial random elongation of the image. As a result, 
both b and its error are large. As s increases, the noise fluctuations decrease and the image 
becomes more circular - both b and its error decrease. At s ~ 10^, the binary nature of 
the star prevents further circularization of the image, and b stabilizes at the correct flux 
ratio. When s ~ 10^, the relative error of b has dropped to about 0.1. From this we can 
conclude that in order to measure the flux ratio with an error of 10%, a total signal-to-noise 
ratio of order 10"^ is required. This is much greater than can be obtained in a single image, 
but could be achieved by coadding many exposures. The exposure time for the individual 
images should be chosen to maximize the signal-to-noise ratio, while remaining in the linear 
region of the detector. As this results in a constraint on the maximum flux in any individual 
pixel, it may be more convenient to work with the peak signal-to-noise ratio Sp, ie. that of 
the central pixel, rather than the total signal-to-noise ratio s. For the photon-noise-limited 
case, the two quantities are related by the square root of the fraction of the total light 
contained within the central pixel 

= (44) 
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Both total and peak signal-to-noise ratios are listed in Table 1. 

The results of the simulation show that the analytic error estimate of Eqn. 43 
is reasonably accurate. For the 1.1-um simulation, this equation predicts that a total 
signal-to-noise ratio s ~ 8000 would be required for a 10% measurement of b which, given 
the approximation in estimating the Variance from Eqn. 28, is in quite good agreement 
with the value s ~ 10000 indicated by the simulation. 

4. DISCUSSION 

We have described a new technique for the detection and measurement of unresolved 
binary systems based on the elongation of their images, in comparison to images of 
refererence stars. The method can only work if variations the shape of the PSF, between 
the binary and reference observations, are smaller than the elongation to be measured. 
The analysis and simulations show what signal-to-noise ratio must be obtained in order to 
reduce the random errors to any desired level. However, one must also consider systematic 
errors that may not be removed by the differential measurement technique. These inchidc 
optical aberrations, guiding errors, and any other effects of this nature. Aberrations can be 
minimized by positioning the reference star and binary at precisely the same location on 
the detector. Guiding errors which are the same for both reference and target observations 
will be removed by the analysis, but any variations will contaminate the signal. Such effects 
should be characterized and understood before observations are attempted. For HST, 
guiding errors are reported to be of order 0.001 arcsec on an individual exposure. If these 
are random, they can be reduced to sufficiently small levels by obtaining many exposures. 

Ideally, it would be best to observe the target and reference stars simultaneously, if 
their separation is sufficiently small that both images can simultaneously fit on the detector. 
In this case asymmetries due to guiding errors should be the same for both images and 
will be removed by the analysis, in order to minimize optical aberrations, the stars should 
be placed on opposite sides of the optical field center, and equidistant from it. To further 
reduce systematic effects, the telescope should be rotated axially by 180 dcg, if possible, to 
interchange the positions of target and reference stars. By exposing for equal times in the 
two configurations, the systematic component of the PSF will then be identical for both 
target and reference images. 

It should be possible to apply this technique to images produced by adaptive optics 
systems on ground-based telescopes. The high-resolution provided by such systems 
facilitates the detection of close binaries, but great care will be needed to minimize 
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variations in the PSF between the binary system and the comparison star. Because the 
performance of adaptive optics systems depends the brightness of the reference star and on 
the atmospheric seeing, the binary system and cahbration star need to be well-matched in 
both magnitude and sky position. We hope to test the feasibility of the method in the near 
future by means of adaptive-optics observations with the Canada-France-Hawaii Telescope. 

I am grateful to Harvey Richer for bringing the problem of the measurement of G77-61 
to my attention, and for several interesting discussions. An anonymous referee provided 
helpful comments and emphasized the importance of matching the comparison and target 
stars. This work was supported by grants from Natural Sciences and Engineering Research 
Council of Canada. 



Table 1. Simulations of the G77-61 Binary System 
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1.10 um 

< b > 




a{b) 


Sp 




1.55 um 

< b > 




a{b) 


10 


4. 


,03 


0.507294 


0. 


,486904 


2. 


,86 


0.472412 


0. 


,494915 


22 


8. 


,68 


0.428661 


0. 


,435401 


6. 


16 


0.460700 


0. 


,476939 


46 


18. 


,71 


0.136120 


0. 


,198945 


13. 


,28 


0.307730 


0. 


,391819 


100 


40. 


,31 


0.049294 


0. 


,065857 


28. 


,61 


0.125452 


0. 


,184577 


215 


86. 


,84 


0.023317 


0. 


,027893 


61. 


,63 


0.054850 


0. 


,067542 


464 


187. 


,10 


0.013118 


0. 


,013616 


132. 


,78 


0.030990 


0. 


,031048 


1000 


403. 


,09 


0.009182 


0. 


,007273 


286. 


,07 


0.022006 


0. 


,016621 


2154 


868. 


,43 


0.008086 


0. 


,003830 


616. 


,30 


0.019957 


0. 


,008461 


4642 


1871. 


,00 


0.007879 


0. 


,001808 


1327. 


,81 


0.019566 


0. 


,003946 


10000 


4030. 


,94 


0.007801 


0. 


,000839 


2860. 


,67 


0.019414 


0. 


,001832 


21544 


8684. 


,26 


0.007766 


0. 


,000389 


6163. 


,02 


0.019347 


0. 


,000850 


46416 


18710, 


.01 


0.007750 


0. 


,000181 


13278. 


,07 


0.019317 


0. 


,000395 


100000 


10309, 


.11 


0.007713 


0, 


.000081 


28606, 


,68 


0.019301 


0, 


,000183 
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Figure Captions 

Fig. 1. — Detecting a simulated binary system. The curves show the estimated flux ratio 
of the binary, and its standard error, as a function of the total signal-to-noise ratio in the 
image. The binary system is assumed to have a true flux ratio of 0.00909 and a separation 
of 0.035 arcsec. The simulation is for a 2.4-meter telescope at a wavelength of 1.10 um. 

Fig. 2. — This plot is similar to Fig. 1, but for a flux ratio of 0.02083 and a wavelength of 
1.55 um. 
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